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A  SIMULATION  STUDY 
OF  RANDOM  CAPS  ON  A  SPHERE 


H.  Solomon  and  C.  Sutton 


Abstract:  This  paper  describes  the  computer  simulation  of  a  coverage  problem  in  geometric  probability, 
that  of  placing  random  caps  on  the  surface  of  a  sphere.  The  simulation  results  were  compared  with 
exact  values  where  known  and  the  differences  were  negligible.  This  suggested  the  use  of  simulation 
results  to  assess  several  approximation  formulas  in  the  literature. 


I 


1.  Introduction 

Consider  a  sphere  of  unit  radius  on  which  are  placed  N  spherical  caps  subtending  a 
half  angle  a  at  the  center  of  the  sphere  (0  <  a  <  x/2).  Assume  that  the  centers  of  these 
caps  are  independently  and  uniformly  distributed  over  the  surface  of  the  sphere.  We  seek  to 
approximate  the  probability  P (N)  that  the  sphere  is  completely  covered.  An  exact  expression 
for  P (N)  is  presently  not  known  except  for  several  special  cases,  see  Solomon  (1978). 

An  approximation  of  P(Af)  can  be  obtained  via  computer  simulation  as  follows.  For  each 
of  n  trials  we  generate  N  random  caps  on  the  surface  of  a  unit  sphere.  We  let  denote 

the  number  of  trials  in  which  the  caps  completely  cover  the  sphere.  Then  pjy,*  =  is 

an  approximate  value  of  P(N).  As  n  increases,  the  approximation  should  approach  the  exact 
value. 

In  order  to  determine  whether  or  not  the  caps  of  a  given  trial  completely  cover  a  sphere, 
let  us  define  a  crossing  as  being  a  point  of  intersection  of  the  circular  boundaries  of  two 
overlapping  caps.  Then  if  there  is  an  uncovered  crossing,  some  area  of  the  sphere  outside 
the  two  overlapping  caps  must  also  be  uncovered.  Conversely,  if  at  least  two  caps  overlap, 
then  the  boundary  of  any  uncovered  area  on  the  sphere  must  contain  an  uncovered  crossing. 
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Thus  we  find  that  the  sphere  is  completely  covered  if  and  only  if  there  are  at  least  two  caps 
which  overlap  and  every  crossing  is  covered.  This  is  the  essential  fact  on  which  the  simulation 


program  is  based. 


Gilbert  (1965)  used  the  preceding  idea  to  perform  a  similar  simulation;  however,  his 
study  was  rather  small  and  the  results  do  not  provide  enough  data  with  which  to  evaluate 
the  accuracy  of  the  various  approximation  formulas  which  have  been  suggested  for  P(JV). 
Prior  results  were  also  achieved  by  Moran  and  Fazekas  de  St.  Groth  (1962)  by  employing  a 
small  physical  simulation.  We  desire  to  assess  the  accuracy  of  the  suggested  approximation 
formulas  and  have  therefore  performed  a  large  computer  simulation.  Some  of  the  main  steps 
incorporated  in  our  simulation  program  are  described  below. 


2.  Simulation 


Random  points  on  the  surface  of  the  sphere,  which  provide  centers  for  the  caps,  are 
obtained  as  follows.  For  each  point  p,-  we  first  generate  2,  and  g,>3  randomly,  each 

according  to  a  normal  distribution  with  mean  0  and  variance  1.  We  then  perform  a  projection 


by  letting 


P,J  "  Hill 


-  ?»\2 

**  “  IMI 

„  9i,i 

P,J ~  IMI 


where 


IMI -M. +*+&)*• 


Due  to  the  symmetry  of  the  spherical  normal  distribution,  the  point  p,-  =  (p.,i,p.,j,Pi,s)  so 
obtained  will  be  a  random  observation  from  a  distribution  which  is  uniformly  distributed  over 
the  surface  of  the  sphere. 


A  cap  intersects  another  cap  if  its  center  falls  within  a  great  circle  distance  2a  from  the 
center  of  the  other  cap.  Letting  d(a,b)  denote  the  straight  line  distance  between  two  points  a 
and  6,  we  then  have  that  the  cap  with  center  a  intersects  the  cap  with  center  b  if  and  only  if 


arcsin(d(a,  b)/2)  <  a. 


Suppose  it  occurs  that  the  cap  with  center  a  =  (01,02,03)  intersects  the  cap  with  cen¬ 
ter  6  =  (61,62,63).  Then  two  crossings  are  formed,  one  on  each  side  of  the  great  circle  arc 
connecting  points  a  and  b.  Denoting  the  midpoint  of  this  arc  by  m,  the  crossings  will  lie  on 
the  great  circle  which  forms  a  perpendicular  intersection  at  the  point  m  with  the  great  circle 
passing  through  a  and  b.  Using  spherical  trigonometry  we  obtain  that  the  great  circle  distance 
between  each  crossing  and  the  point  m  is  given  by 


arccos  < 


.cos  arcsin 


so  1  I  cos  a  | 

— t,  rr~  t  =  arccos  < - r  >. 

in^]/  l[,-Ci±*l]!j 


Denote  this  distance  by  7,  and  call  the  two  crossings  formed  z  and  y. 

We  can  determine  the  Cartesian  coordinates  (zi,z2,zs)  of  z  and  (yi,V2,Vi)  of  y  in  the 
following  manner  Set  v  =  (vi,«2,«s)  equal  to  (61  -  ai,&2  —  02,63  -  03).  Let  ||v||  denote 
the  length  of  the  vector  v.  Let  6  be  the  angle  between  the  positive  x-axis  and  the  vector 
(vi,V2,0),  measured  in  the  direction  from  the  positive  x-axis  towards  the  positive  y-axis.  Let 
4>  =  arcsin( V3/I |v| |).  Consider  the  following  series  of  rotation  of  axes.  First  rotate  about  the 
s-axis  an  angle  6  in  the  direction  of  (wi,»2>0)-  Then  rotate  about  the  new  y-axis  an  angle 
<f>  towards  the  positive  z-axis.  Next  we  rotate  about  the  new  x-axis  an  angle  7  towards  the 
positive  y-axis.  Now  if  we  reverse  the  direction  of  the  first  two  rotations,  and  rotate  about  the 
y-axis  by  <j>  and  then  about  the  z-axis  by  0,  the  new  coordinates  of  m  will  be  the  coordinates 
of  the  crossing  z  in  the  original  system.  This  sequence  of  rotations  of  axes  is  equivalent  to 
rotating  the  sphere  through  an  angle  7  about  an  axis  that  passes  through  the  center  of  the 
sphere  and  which  is  parallel  to  the  vector  v.  The  coordinates  of  y  are  found  in  the  same  way, 
except  that  we  rotate  about  the  x-axis  by  an  angle  -7. 

A  crossing  formed  by  two  caps  will  be  covered  if  the  great  circle  distance  between  the 
crossing  and  the  center  of  any  of  the  N- 2  other  caps  is  less  than  a.  Equivalently,  we  can  check 
whether  or  not  the  straight  line  distance  is  less  than  2sin(a/2). 

Simulations  were  done  using  a  =  v/2  and  N  =  4, 5, . . . ,  13.  This  value  of  a  was  selected 
since  exact  values  for  P (N)  have  been  determined  in  this  case  and  hence  can  serve  as  anchors 
in  our  simulation  study.  J.  G.  Wendel  (1962)  showed  that  if  N  points  are  scattered  at  random 
on  the  surface  of  the  unit  sphere  in  n-space,  the  probability  that  all  the  points  lie  on  some 


hemisphere  is  given  by 


2  (';*)■ 

i=0  '  ' 

This  expression  yields,  as  a  special  case,  the  result  that  for  spherical  caps  of  angular  radius 
a  =  jr/2 


»-iV+l 


(1) 


P (N)  =  1  -  2 ~n(N2  -N  +  2). 


For  each  value  of  N  20, (MX)  trials  were  accomplished.  The  number  of  covered  spheres  out 
of  20,000  trials,  tojv^o ooo.  i*  binomially  distributed  with  parameters  n  =  20,000  and  p  =  P (N). 
Hence  the  standard  deviation  of  pjv, 20000  ( =  **1x20000/ 2 0000 )  is  given  by 

SDN,mooo  -  V - 20OOO - • 

The  following  results  give,  for  each  value  of  N,  the  exact  value  of  P(N)  for  a  =  jt/2  ,  our 
approximation  pn, 20000,  and  the  value  of  1x20000  =  (pnjoooo  ~  P[N))/SDff<2oooo- 


N 

P(JV) 

PX20000 

tn, 20000 

4 

0.12500 

0.12550 

0.21 

5 

0.31250 

0.31555 

0.93 

6 

0.50000 

0.50030 

0.08 

7 

0.65625 

0.66055 

1.28 

8 

0.77344 

0.77145 

-0.67 

0 

0.85547 

0.85720 

0.70 

10 

0.91016 

0.91100 

0.42 

11 

0.94531 

0.94430 

-0.63 

12 

0.96729 

0.96825 

0.76 

13 

0.98071 

0.98090 

0.20 

These  results  serve  to  validify  our  simulation  program. 


S.  Moran  and  Fazekas  de  St.  Groth  result* 

We  desire  to  investigate  the  various  formulas  that  have  been  developed  to  approximate 
P(1V).  First  the  formulas  will  be  compared  for  the  case  of  a  =  53.43*.  This  value  of  a  is  chosen 
since  it  has  been  considered  by  previous  investigators  (because  of  a  biological  application).  To 
serve  as  a  basis  for  the  comparisons  we  estimated  P(1V)  using  our  simulation  program  for 
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N  =  10,  IS,  20, . . . ,  60 (  and  a  =  53.43°).  Each  estimate  is  based  on  20,000  trials.  We  note  that 
even  for  20,000  trials  the  standard  deviation  of  the  estimator  pw, 20000  can  be  as  large  as  0.0035 
for  some  values  of  N.  Hence  the  displayed  simulation  results  in  this  paper  may  indicate  more 
accuracy  than  is  actually  present. 

Moran  and  Fazekas  de  St.Groth  (1962)  used  a  method  of  moments  approach  to  derive 
two  approximation  formulas  for  P (N).  Under  the  assumption  that  the  uncovered  region  of  the 
sphere  consists  of  a  single  region  (which  will  always  be  the  case  when  a  =  jr/2  and  should  be 
nearly  true  for  a  close  to  x/2)  they  developed 

(2)  P(N)~l-i»2[£(Y)]2/S(r*), 

where  Y  denotes  the  proportion  of  the  surface  not  covered.  It  is  found  that 

B(Y)  =  cosJAr(a/2) 

and  the  calculation  of  E(Y2)  may  be  accomplished  via  numerical  integration  using 

(3)  E(yJ)=  *jf[l"  ^/(*)]WsM«f* 
and 

(4)  1  _  Ml  =  H1  "  »  arcco8(*^^i)] cosa  +  i  arccos(^ ?£*)  (°  £  *  £  2“) 

(cosa  (2a  <  ^  <  jr). 

Under  the  hypothesis  that  the  number  of  disjoint  uncovered  regions  has  a  Poisson  distri¬ 
bution  and  that  the  areas  of  the  individual  uncovered  regions  are  distributed  independently, 
Moran  and  Fazekas  de  St.  Groth  derived 

(5)  P(JV)  =  exp  (-^2{E(y2)[E(y)]-J  -  1}-1). 

They  expected  that  the  actual  value  of  P (N)  should  be  between  the  values  given  by  (2)  and 

(5). 

In  the  table  below  we  give  our  simulation  results  for  P(AT)  for  a  =  53.43°  based  on  20,000 
trials  for  each  value  of  N,  as  well  as  the  values  obtained  from  both  of  the  formulas  of  Moran 
and  Fazekas  de  St.  Groth. 
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N  Simulation  (2),(3)&(4)  (5),(3)&(4) 


10 

0.0003 

-2.6119 

0.0000 

15 

0.0465 

-1.3333 

0.0120 

20 

0.2752 

-0.2882 

0.1750 

25 

0.5763 

0.3643 

0.4821 

30 

0.7846 

0.7085 

0.7336 

35 

0.8903 

0.8738 

0.8785 

40 

0.9556 

0.9474 

0.9482 

45 

0.9837 

0.9787 

0.9788 

50 

0.9929 

0.9915 

0.9916 

55 

0.9974 

0.9967 

0.9967 

60 

0.9988 

0.9987 

0.9987 

It  appears  that  both  approximation  formulas  converge  to  a  common  value  as  N  increases; 
however,  they  do  not  sandwich  the  simulation  value  as  was  anticipated  by  Moran  and  Fazekas 
de  St.  Groth.  Rather,  both  formulas  tend  to  underestimate  P(AT). 


In  their  paper,  Moran  and  Fazekas  de  St.  Groth  did  not  use  (3)  and  (4)  in  their  calcula¬ 
tions  of  P(JV).  Instead  they  use  a  saddlepoint  method  to  arrive  at 

(6)  = 

as  an  approximation  of  E(Y2).  It  is  interesting  to  note  that  making  use  of  this  substitution 
leads  to  values  of  P (N)  which  better  approximate  these  simulation  values.  A  comparison  is 
presented  in  the  table  below,  where  for  each  case  it  can  be  seen  that  the  approximation  (6) 
yields  a  better  estimate. 


N 

Simulation 

Value  from 

Value  from 

Value  from 

Value  from 

value 

(2),(3)&(4) 

(2)  and  (6) 

(5),(3)&(4)  (5)  and  (6) 

10 

0.0003 

•2.6119 

-2.6814 

0.0000 

0.0000 

15 

0.0465 

-1.3333 

-1.2616 

0.0120 

0.0154 

20 

0.2752 

•0.2882 

•0.2163 

0.1750 

0.1991 

25 

0.5763 

0.3643 

0.4050 

0.4821 

0.5084 

30 

0.7846 

0.7085 

0.7279 

0.7336 

0.7498 

35 

0.8993 

0.8787 

0.8815 

0.8785 

0.8857 

40 

0.9556 

0.9474 

0.9503 

0.9482 

0.9511 

45 

0.9837 

0.9787 

0.9798 

0.9788 

0.9799 

50 

0.9929 

0.9915 

0.9920 

0.9916 

0.9920 

55 

0.9974 

0.9967 

0.9969 

0.9967 

0.9969 

60 

0.9988 

0.9987 

0.9988 

0.9987 

0.9988 
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For  a  =  53.43°  and  N  =  10, 15,  ...,60  we  found  that  the  combination  of  (5)  and  (6)  constituted 
the  best  approximation  for  P (N)  produced  by  Moran  and  Fazekas  de  St.  Groth. 

4.  Gilbert  developement 

E.  N.  Gilbert  (1965)  developed  upper  and  lower  bounds  for  P(./V).  Let  X  denote  the 
probability  that  any  specified  point  on  the  sphere  will  be  covered  by  a  randomly  placed  cap, 
so  that 

X  =  ^(1  -  cos  a)  =  sin2(a/2). 

m 

The  lower  bound  for  P (N)  is  given  by 

(7)  1-±N(N-1)X(1-X)n-1. 

Gilbert  also  gives 

(8)  l-(l-A)* 

as  a  general  upper  bound  for  P (N).  This  is  a  weak  upper  bound.  For  the  specific  case  of 
a  =  53.43°  a  better  upper  bound  can  be  found.  The  dominant  terms  of  this  closer  upper 
bound  are 

(9)  1  -  6(1  -  X)N. 

Note  that  (8)  is  just  the  probability  that  some  fixed  point  V  on  the  sphere  is  covered,  and 
notice  that  the  probability  that  V  is  covered  exceeds  the  probability  that  the  entire  sphere 
is  covered.  To  arrive  at  (9),  consider  points  Vj,V2,  ...,V*  which  are  points  of  intersection  of 
the  sphere  and  of  a  regular  octahedron  inscribed  within  the  sphere.  The  method  of  inclusion 
and  exclusion  can  be  used  to  find  an  exact  expression  for  the  probability  that  the  set  of  six 
vertices  is  covered  by  N  .  This  probability  will  be  an  upper  bound  for  P (N).  For  large  N,  (9) 
gives  the  dominant  terms  of  this  expression.  It  can  be  seen  from  the  table  below  that  while 
the  upper  and  lower  bounds  do  indeed  bracket  the  simulation  values,  neither  provides  a  very 
good  estimate  of  P{N). 
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‘/V? 


HV.‘ y 


N 


lower  bound  (7) 


simulation 


upper  bound  (9) 


10 

-2.1790 

0.0003 

0.3725 

15 

•1.3989 

0.0465 

0.7970 

20 

•0.4039 

0.2752 

0.9344 

25 

0.2831 

0.5763 

0.9788 

30 

0.6638 

0.7846 

0.9931 

35 

0.8513 

0.8993 

0.9978 

40 

0.9370 

0.9556 

0.9993 

45 

0.9741 

0.9837 

0.9998 

50 

0.9896 

0.9929 

0.9999 

55 

0.9959 

0.9974 

1.0000 

60 

0.9984 

0.9988 

1.0000 

5.  Comparison  of  approximations  for  a  =  53.43° 

R.  E.  Miles  (1969)  provides  yet  another  way  of  approximating  P (N).  Theorem  3  of  his 
paper  yields,  as  a  special  case,  the  approximation 

(10)  1  -  2~n N(N  -  1)  sin2  o(l  +  cos a)^-2 


for  P(1V).  The  table  below  compares  values  for  a  =  53.43°  obtained  from  (10)  with  our 
simulation  results.  The  best  estimates  of  Moran  and  Fazekas  de  St.  Groth,  using  (5)  and  (6), 
and  the  best  estimates  of  Gilbert,  using  his  lower  bound  (7),  are  included  for  comparison. 


N 

simulation 

result 

Miles 

Moran  and 

Fazekas  de  St.  Grotb 

Gilbert 

10 

0.0003 

-1.3842 

0.0000 

-2.1790 

15 

0.0465 

-0.7992 

0.0154 

-1.3989 

20 

0.2752 

-0.0529 

0.1991 

-0.4039 

25 

0.5763 

0.4623 

0.5084 

0.2831 

30 

0.7846 

0.7479 

0.7498 

0.6638 

35 

0.8993 

0.8885 

0.8857 

0.8513 

40 

0.9556 

0.9527 

0.9511 

0.9370 

45 

0.9837 

0.9806 

0.9799 

0.9741 

50 

0.9929 

0.9922 

0.9920 

0.9896 

55 

0.9974 

0.9970 

0.9969 

0.9959 

60 

0.9988 

0.9988 

0.9988 

0.9984 

We  can  see  that  although  Moran  and  Fazekas  de  St.  Groth’s  candidate  does  the  best  for  the 
smaller  values  of  N,  Miles’  formula  seems  to  be  the  most  accurate  for  the  higher  values  of  N. 
We  note  that  each  of  the  formulas  appear  to  underestimate  P{N). 
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0.  Comparison  of  approximations  for  a  =  r/2 


In  order  to  provide  another  basis  of  comparison  for  the  various  approximation  schemes, 
we  have  also  calculated  each  method’s  estimate  of  P(Af)  for  a  —  v/2  and  N  =  4,5,6,...,  25. 
We  can  judge  the  accuracy  of  the  estimates  by  comparing  them  with  the  exact  value  of  P(W) 
obtained  from  (1). 

Let  us  first  examine  the  four  candidates  from  the  paper  of  Moran  and  Fasekas  de  St. 
Groth.  Recall  that  they  derived  (2)  and  (5)  from  different  hypotheses,  and  that  for  each 
formula  we  can  either  use  the  exact  value  of  E(Y2)  obtained  by  numerical  integration  or  we 
can  make  use  of  the  approximation  (6).  For  N  =  5, 10, 15, 20,25  (and  a  =  w/2)  the  estimates 
of  P(7V)  from  these  four  methods,  as  well  as  the  known  values  of  P (N),  are  displayed  in  the 
table  below. 


N 

W) 

Value  from 
(2),(3)&(4) 

Value  from 
(2)  and  (6) 

Value  from 

(5),(3)&(4) 

Value  from 
(5)  and  (6) 

5 

0.3125 

-0.5009 

-0.0897 

0.1157 

0.2470 

10 

0.9102 

0.8640 

0.8927 

0.8695 

0.8961 

15 

0.9935 

0.9915 

0.9928 

0.9915 

0.9928 

20 

0.9996 

0.9996 

0.9996 

0.9996 

0.9996 

25 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

One  can  see  that  for  a  =  x/2,  the  best  estimate  for  ?(N)  comes  from  using  (5)  and  (6).  This 
was  also  the  case  for  a  =  53.43°.  Also  note  that  each  approximation  formula  underestimates 
P(7V),  as  was  the  case  previously.  Because  there  is  only  one  uncovered  region  when  a  =  v/2,  it 
is  somewhat  surprising  that  (5)  approximates  ?{N)  better  than  does  (2),  since  (2)  was  derived 
under  the  hypothesis  that  there  was  only  one  vacant  region. 

For  a  =  x/2  we  can  again  use  (7)  to  obtain  a  lower  bound  for  P (N),  and  we  can  use  (8) 
to  obtain  an  upper  bound  for  P (N).  These  bounds,  both  by  Gilbert,  and  the  exact  value  of 
P(iV)  are  given  below  for  various  values  of  N. 


N 

lower  bound  (7) 

known  value  (1) 

upper  bound  (8) 

5 

0.1667 

0.3125 

0.9688 

10 

0.8828 

0.9102 

0.9990 

15 

0.9915 

0.9935 

1.0000 

20 

0.9995 

0.9996 

1.0000 

25 

1.0000 

1.0000 

1.0000 
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In  the  table  below  we  display  the  approximated  values  of  P (N)  derived  from  the  work 
of  Miles,  equation  (10)  Also  shown  are  the  best  approximations  of  Moran  and  Fazekas  de 
St  Groth,  (5)  used  with  (6),  and  Gilbert’s  lower  bound  (7).  The  known  values  of  P(N)  from 
equation  (1)  are  also  given. 


N 

known 

value 

Miles 

Moran  and 

Fazekas  de  St.  Groth 

Gilbert 

4 

0  1250 

0.2500 

0.0903 

0.0000 

6 

0  5000 

0.5313 

0.4324 

0.3750 

8 

0.7734 

0.7813 

0.73G0 

0.7083 

10 

09102 

0.9121 

0.8961 

0.8828 

12 

09673 

0.9678 

0.9629 

0.9570 

14 

09888 

0.9889 

0.9875 

0.9852 

16 

0.9963 

0.9963 

0.9959 

0.9951 

18 

0  9988 

0.9988 

0.9987 

0.9984 

20 

0.9996 

0.9996 

0.9996 

0.9995 

22 

0  9999 

0.9999 

0.9999 

0.9999 

24 

1.0000 

1.0000 

1.0000 

1.0000 

It  can  be  seen  that  the  formula  of  Miles  provides  the  most  accurate  estimates  of  P (N),  and 
also  that  Moran  and  Fazekas  de  St.  Groth’s  candidate  is  better  than  Gilbert’s.  Note  that 
Miles’  formula  overestimates  P(N)  for  this  value  of  a,  whereas  before,  with  a  =  53.43°,  it 
underestimated  P(N)  It  is  interesting  to  note  that  for  a  =  x/2  Miles’  formula  reduces  to 

P(N)  n  1  -  2 ~n(N7  -  N) 
which  is  very  similar  to  the  exact  result, 

P(!V)  =  1  -  2 ~S(N7  -  N+  2). 
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